﻿#pragma once
#include <random>
#include <algorithm>
#include <vector>
#include <iostream>
#include <iomanip>
#include "Vec3.h"
#include "util.h"


//from ray tracing:The rest of your life
struct sample {
	double x;
	double p_x;
};

bool compare_by_x(const sample& a, const sample& b) {
	return a.x < b.x;
}


static void half_area() {
	const unsigned int N = 10000;
	sample samples[N];
	double sum = 0.0;

	// Iterate through all of our samples.

	for (unsigned int i = 0; i < N; i++) {
		// Get the area under the curve.
		auto x = random_double(0, 2 * pi);
		auto sin_x = std::sin(x);
		auto p_x = exp(-x / (2 * pi)) * sin_x * sin_x;
		sum += p_x;

		// Store this sample.
		sample this_sample = { x, p_x };
		samples[i] = this_sample;
	}

	// Sort the samples by x.
	std::sort(std::begin(samples), std::end(samples), compare_by_x);

	// Find out the sample at which we have half of our area.
	double half_sum = sum / 2.0;
	double halfway_point = 0.0;
	double accum = 0.0;
	for (unsigned int i = 0; i < N; i++) {
		accum += samples[i].p_x;
		if (accum >= half_sum) {
			halfway_point = samples[i].x;
			break;
		}
	}

	std::cout << std::fixed << std::setprecision(12);
	std::cout << "Average = " << sum / N << '\n';
	std::cout << "Area under curve = " << 2 * pi * sum / N << '\n';
	std::cout << "Halfway = " << halfway_point << '\n';
}

/*

 CDF: Cumulative Distribution Function
 ICD: inverse of our CDF

 solution the Integral of $int_{0}^{2}x^2 dx $
 uniform sample between [0,2],then mean that we use a PDF = 1/2 over the rangle [0-2],
 which mean the CDF is P(x) = x/2,so the ICD(d) = 2d.
*/

void integral_power_of_x_with_uniform_sample() {

	//this equal the f(d)
	auto icd = [](double d) {
		return 2.0 * d;
	};

	auto pdf = [](double x) {
		return 0.5;
	};

	//integral  $int_{0}^{x}x^2 dx 
	auto f = [](double x) {
		return x * x;
	};

	int N = 1000000;
	auto sum = 0.0;

	for (int i = 0; i < N; i++) {

		auto d = icd(random_double());

		sum += f(d) / pdf(d);
	}
	std::cout << std::fixed << std::setprecision(12);
	std::cout << "I = " << (sum / N) << '\n';
}

/*

We need to account for the nonuniformity of the PDF of x
we need to down-weight where we sample more frequently, and to up-weight where we sample less frequently.
For our new nonuniform random number generator, the PDF defines how much or how little we sample a specific portion.
So the weighting function should be proportional to 1/pdf
. In fact it is exactly 1/pdf

We can try to solve for the integral using the linear PDF, p(r)=r/2.
To do that, all we need to do is replace the functions $ICD(d) = \sqrt(d4)$
 and p(x)=x/2  任何一个pdf都可以，只是这里使用了这个不是uniformity 的pdf
*/

void integral_power_of_x_with_nonuniformity_sample() {

	auto icd = [](double d) {
		return std::sqrt(4 * d);
	};

	auto pdf = [](double x) {
		return x * 0.5;
	};

	//integral  $int_{0}^{x}x^2 dx 
	auto f = [](double x) {
		return x * x;
	};

	int N = 1000000;
	auto sum = 0.0;

	for (int i = 0; i < N; i++) {
		auto z = random_double();
		if (z == .0) continue;
		auto d = icd(z);
		sum += f(d) / pdf(d);
	}
	std::cout << std::fixed << std::setprecision(12);
	std::cout << "I = " << (sum / N) << '\n';
}

/*
	we should just try to make the PDF match the integrand by turning the PDF into a quadratic function:

 

In all of the examples given, we always converged to the correct answer of 8/3
. We got the same answer when we used both a uniform PDF and the “correct” PDF (that is, ICD(d)=8d13
). While they both converged to the same answer, the uniform PDF took much longer.
After all, we only needed a single sample from the PDF that perfectly matched the integral.
This should make sense, as we were choosing to sample the important parts of the distribution more often,
whereas the uniform PDF just sampled the whole distribution equally, without taking importance into account.


*/

void integral_power_of_x_with_equal_function_sample() {

	auto icd = [](double d) {
		return 8.0 * std::pow(d, 1.0 / 3.0);
	};

	auto pdf = [](double x) {
		return (3.0 / 8.0) * x * x;
	};

	//integral  $int_{0}^{x}x^2 dx 
	auto f = [](double x) {
		return x * x;
	};

	int N = 1;
	auto sum = 0.0;

	for (int i = 0; i < N; i++) {
		auto z = random_double();
		if (z == 0.0)  // Ignore zero to avoid NaNs
			continue;

		auto d = icd(z);
		sum += f(d) / pdf(d);
	}
	std::cout << std::fixed << std::setprecision(12);
	std::cout << "I = " << sum / N << '\n';
}



/*
we integration : f(\theta,\phi) = cos^2(\theta)  in unit sphere

Using Monte Carlo integration, we should just be able to sample cos^2(θ)/p(r)
What is the PDF of these uniform samples  it is 1/area of the sphere, which is 1/(4π)
*/

void integral_over_sphere()
{
	auto pdf = [](const Vec3& d) {
		return 1 / (4 * pi);
	};

	auto  f = [](const Vec3& d) {
		auto cosine_squared = d.z() * d.z();
		return cosine_squared;
	};

	int N = 1000000;
	auto sum = 0.0;
	for (int i = 0; i < N; i++) {
		Vec3 d = random_unit_vector();
		auto f_d = f(d);
		sum += f_d / pdf(d);
	}
	std::cout << std::fixed << std::setprecision(12);
	std::cout << "I = " << sum / N << '\n';
}


/////////////////////////////////////////////////////////////////////////////
//from ray tracing:The rest of your life
Vec3 generate_random_on_sphere() 
{
	for (int i = 0; i < 200; i++) {
		auto r1 = random_double();
		auto r2 = random_double();
		auto x = std::cos(2 * pi * r1) * 2 * std::sqrt(r2 * (1 - r2));
		auto y = std::sin(2 * pi * r1) * 2 * std::sqrt(r2 * (1 - r2));
		auto z = 1 - 2 * r2;

		return Vec3(x, y, z);
	}

	return Vec3(.0, .0, .0);
}
//from ray tracing:The rest of your life
Vec3 generate_random_cosine_direction() {
	auto r1 = random_double();
	auto r2 = random_double();

	auto phi = 2 * pi * r1;
	auto x = std::cos(phi) * std::sqrt(r2);
	auto y = std::sin(phi) * std::sqrt(r2);
	auto z = std::sqrt(1 - r2);

	return Vec3(x, y, z);
}


/*
*
*Integrate f(θ) = cos(θ)^3 over  Hemisphere
* Now for integration with importance sampling. p(ω)=1/2π
, so we average f()/p()=cos(θ)^3/1/2π
*/

void  integration_over_Hemisphere() {

	auto pdf = []() {
		return 1.0 / (2.0 * pi);
	};

	auto  f = [](double r2) {
		auto z = 1 - r2;
		double cos_theta = z;
		return cos_theta * cos_theta * cos_theta;
	};

	int N = 1000000;
	auto sum = 0.0;
	for (int i = 0; i < N; i++) {
		auto r2 = random_double();
		sum += f(r2) / pdf();
	}

	std::cout << std::fixed << std::setprecision(12);
	std::cout << "PI/2 = " << pi / 2.0 << '\n';
	std::cout << "Estimate = " << sum / N << '\n';
}


void integration_over_Hemisphere_with_consine() {

	auto  f = [](const Vec3& d) {
		auto cos_theta = d.z();
		return cos_theta * cos_theta * cos_theta;
	};

	auto pdf = [](const Vec3& d) {
		return d.z() / pi;
	};

	int N = 1000000;

	auto sum = 0.0;
	for (int i = 0; i < N; i++) {
		Vec3 d = generate_random_cosine_direction();
		sum += f(d) / pdf(d);
	}

	std::cout << std::fixed << std::setprecision(12);
	std::cout << "PI/2 = " << pi / 2.0 << '\n';
	std::cout << "Estimate = " << sum / N << '\n';
}


/*
*
If the light does scatter, it will have a directional distribution that we can describe as a PDF over solid angle.

I will refer to this as its scattering PDF: pScatter()
- The scattering PDF will vary with outgoing direction: pScatter(ωo)
- The scattering PDF can also vary with incident direction: pScatter(ωi,ωo)
- The scattering PDF can vary with the wavelength of the light: pScatter(ωi,ωo,λ)

The albedo of an object can also depend on these quantities: A(x,ωi,ωo,λ)

The color of a surface is found by integrating these terms over the unit hemisphere by the incident direction:
Coloro(x,ωo,λ)=∫ωiA(x,ωi,ωo,λ)⋅pScatter(x,ωi,ωo,λ)⋅Colori(x,ωi,λ)

We've added a Colori term. The scattering PDF and the albedo at the surface of an object are acting as filters to
 the light that is shining on that point. So we need to solve for the light that is shining on that point.


 For a Lambertian surface we already implicitly implemented this formula for the special case where pScatter(…) is a cosine density.
 The pScatter(…) of a Lambertian surface is proportional to cos(θo), where θo is the angle relative to the surface normal. Let's solve for C
 once more:
	pScatter(x,ωi,ωo,λ)=C⋅cos(θo)

	then $\int_{0}^{2\pi}\int_{0}^{\pi} C * cos(θo)dA = 1$  ====> C = 1/pi
	so pScatter(…) = cos(θo) / pi  (we can see this in sample_f method of lambertian )

 */

 /*
  Import sample

 Our goal over the next several chapters is to instrument our program to send a bunch of extra rays toward light sources
 so that our picture is less noisy.

 Let’s assume we can send a bunch of rays toward the light source using a PDF pLight(ωo)
 . Let’s also assume we have a PDF related to pScatter, and let’s call that pSurface(ωo)
 . A great thing about PDFs is that you can just use linear mixtures of them to form mixture densities that are also PDFs. For example, the simplest would be:


 p(ωo)=1/2 pSurface(ωo)+ 1/2 pLight(ωo)
 As long as the weights are positive and add up to one, any such mixture of PDFs is a PDF.
 Remember, we can use any PDF: all PDFs eventually converge to the correct answer.
 So, the game is to figure out how to make the PDF larger where the product

 pScatter(x,ωi,ωo)⋅Colori(x,ωi)
 is largest. For diffuse surfaces, this is mainly a matter of guessing where Colori(x,ωi)
  is largest. Which is equivalent to guessing where the most light is coming from.

 For a mirror, pScatter()
  is huge only near one direction, so pScatter()
  matters a lot more. In fact, most renderers just make mirrors a special case, and make the pScatter()/p()
  implicit — our code currently does that.


 */

 /// <summary>
 /// We have used a PDF related to cos(θ)
 //  and a PDF related to sampling the light.We would like a PDF that combines these.
 /// </summary>
 /// <param name="pWorld"></param>
 /*
 color Camera::ray_color(const Ray& r, int depth,const Hittable& world) const
 {
	 if (depth < 0) return color(0, 0, 0);

	 HitRecord rec;
	 if (!world.hit(r, Interval(0.001, std::numeric_limits<double>::infinity()), rec)) {
		 return background;
	 }

	 Ray scattered;
	 color attenuation;
	 double pdf;

	 color color_from_emission = rec.material->emitted(r,rec,rec.u, rec.v, rec.p);

	 if (!rec.material->scatter(r, rec, attenuation, scattered, pdf))
		 return  color_from_emission;

	 if (!scattering_pdf_of_light(rec, pdf, scattered, r.time()))
		 return color_from_emission;

	 double scattering_pdf = rec.material->scattering_pdf(r, rec, scattered);

	 color color_from_scatter =
		 (attenuation * scattering_pdf * ray_color(scattered, depth - 1, world)) / pdf;

	 //color color_from_scatter = attenuation * ray_color(scattered, depth - 1, world);

	 return color_from_emission + color_from_scatter;


  }
 */

void testx() {
	// view plane
	//half_area();
	integral_power_of_x_with_uniform_sample();

	integral_power_of_x_with_nonuniformity_sample();

	integral_power_of_x_with_equal_function_sample();

	integral_over_sphere();

	integration_over_Hemisphere();

	integration_over_Hemisphere_with_consine();
}


